Тренд

library("TSA")
## Loading required package: leaps
## Loading required package: locfit
## locfit 1.5-9.1    2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.
## Loading required package: tseries
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar

Понятие тренда имеет довольно субъективный характер. В частности, в процессе случайного блуждания можно сказать, что присутствует тренд

n<-200
e <- rnorm(n)
x <- diffinv(e,lag = 1,differences = 1,0)
matplot(x,type = "l",main= "Simulated random walk",lwd = 3,col = "blue",xlab ="time",ylab = "random walk")

Однако известно, что \(E[x_t]=0\) для любого \(t\). Некоторые авторы их называют такие случайные тренды “стохастическими трендами”.

Постоянное среднее

Начнем с простейшей модели

\(y_t= x_t+\mu: t=i,...,n\)

где \(E[x_t]=0\) для всех \(t\)

Чаше всего в качестве оценки \(\mu\) используют обычное среднее

\(\mu=\widehat{y}= \frac{\sum_{i=1}^n}{n}y_i\)

Пусть \(y_t\) или что эквиалентно \(x_t\) стационарный процесс с автокорреляционной функцией \(\rho_k\). Тогда дисперсия оценки \(D[\widehat{y_t}]\) равна

\(D[\widehat{y_t}]=\frac{c_0}{n}\sum_{k=-n+1}^{n+1}(1-\frac{|k|}{n})\rho_k\)

Для нестационарного процесса случайного блуждания \(x_t\)

\(D[\widehat{y_t}]=\frac{1}{n^2}D[\sum_{i=1}^nY_i]=\) \(\frac{1}{n^2}D[\sum_{i=1}^n\sum_{j=1}^i\epsilon_j]=\) \(\frac{1}{n^2}D[\epsilon_1+2\epsilon_2+3\epsilon_3+...+n\epsilon_n]=\) \(\frac{1}{n^2}\sigma_{\epsilon}^2\sum_{k=1}^nk^2=\) \(\sigma_{\epsilon}^2(2n+1)\frac{n+1}{6n}\)

Для произвольного стационарного процесса с автокорреляционной функцией \(\rho_k\) такой, что

\(\sum_{k=1}^n|\rho_k|<\infty\)

При \(n->\infty\) справедливо следующее приближение

\(D[\widehat{y_t}]\approx\frac{c_0}{n}\sum_{k=-\infty}^\infty\rho_k\) (1)

В качестве примера пусть у процесса \(x_t\) автокорреляционная функция \(\rho_k=\phi^k\) тогда

\(D[\widehat{y}_t]\approx\frac{(1+\phi)*c_0}{n(1-\phi)}\)

Пример такого процесса

\(x_t=\phi x_{t-1}+\epsilon_t\) при \(|\phi|<1\)

Это процесс называется процессом авторегрессии первого порядка.

phi <- 0.9
ar1 <- vector()
ar1 <- c(ar1,e[1])
for (i in 2:n)
{
  ar1 <- c(ar1,phi*ar1[i-1]+e[i])
}   
matplot(ar1,main = "x(t)=phi*x(t-1)+e(t)", col="blue",type="l",lwd=3)

acf(ar1,type= "correlation")

Простая линейная регрессия

Пусть \(\mu_t\) неслучайный тренд вида

\(\mu_t=\beta_0+\beta_1t: t=1,2,...,n\)

Минимизируюя сумму квадратов ошибок

\(RSS(\beta_0,\beta_1) = \sum_{t=1}^n(y_t-(\beta_0+\beta_1t))^2\)

найдем оценки параметров \(\beta_0\) и \(\beta_1\). Вычислим частные производный и, приравнивая их нулю, получим

\(\frac{\partial RSS(\beta_0,\beta_1)}{\partial{\beta_0}}= 0\)

\(\frac{\partial RSS(\beta_0,\beta_1)}{\partial{\beta_1}}= 0\)

получим

\(\widehat{\beta}_1=\frac{\sum_{t=1}^n(y_t-\overline y)(t-\overline t)}{\sum_{t=1}^n(t-\overline t)^2}\)

\(\beta_0=\overline{y} -\beta_1t\)

здесь \(\overline{t}=\frac{n+1}{2}\)

Пример для случайного блуждания

model1<- lm(x ~time(x))
summary(model1)
## 
## Call:
## lm(formula = x ~ time(x))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.041  -2.301   1.514   4.205   8.890 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.172910   0.789024  -6.556 4.65e-10 ***
## time(x)     -0.058242   0.006774  -8.598 2.40e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.572 on 199 degrees of freedom
## Multiple R-squared:  0.2709, Adjusted R-squared:  0.2672 
## F-statistic: 73.92 on 1 and 199 DF,  p-value: 2.396e-15

Позднее в этой лекции дадим описание и назначение всех статистик в этом выводе

plot(x , col="blue" , type ="l", lwd=3)
abline(model1)

Циклический тренд

Предположим модель будет

\(y_t=\mu_t+x_t\)

где \(E[x_t]=0\),а \(\mu_t\) сезонные данные т.е. 12 констант \(\beta_1,\beta_2,...,\beta_{12}\) Это модель иногда называют моделью сезонного среднего

library("TSA")
data(tempdub)
plot(tempdub)

Сначала подгоним модель без константы

month.<- season(tempdub)
model2 <- lm(tempdub~month.-1)
summary(model2)
## 
## Call:
## lm(formula = tempdub ~ month. - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2750 -2.2479  0.1125  1.8896  9.8250 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## month.January     16.608      0.987   16.83   <2e-16 ***
## month.February    20.650      0.987   20.92   <2e-16 ***
## month.March       32.475      0.987   32.90   <2e-16 ***
## month.April       46.525      0.987   47.14   <2e-16 ***
## month.May         58.092      0.987   58.86   <2e-16 ***
## month.June        67.500      0.987   68.39   <2e-16 ***
## month.July        71.717      0.987   72.66   <2e-16 ***
## month.August      69.333      0.987   70.25   <2e-16 ***
## month.September   61.025      0.987   61.83   <2e-16 ***
## month.October     50.975      0.987   51.65   <2e-16 ***
## month.November    36.650      0.987   37.13   <2e-16 ***
## month.December    23.642      0.987   23.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.419 on 132 degrees of freedom
## Multiple R-squared:  0.9957, Adjusted R-squared:  0.9953 
## F-statistic:  2569 on 12 and 132 DF,  p-value: < 2.2e-16

Посмотрим как изменится модель если включить оценку константы

month.<- season(tempdub)
model3 <- lm(tempdub~month.)
summary(model3)
## 
## Call:
## lm(formula = tempdub ~ month.)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2750 -2.2479  0.1125  1.8896  9.8250 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       16.608      0.987  16.828  < 2e-16 ***
## month.February     4.042      1.396   2.896  0.00443 ** 
## month.March       15.867      1.396  11.368  < 2e-16 ***
## month.April       29.917      1.396  21.434  < 2e-16 ***
## month.May         41.483      1.396  29.721  < 2e-16 ***
## month.June        50.892      1.396  36.461  < 2e-16 ***
## month.July        55.108      1.396  39.482  < 2e-16 ***
## month.August      52.725      1.396  37.775  < 2e-16 ***
## month.September   44.417      1.396  31.822  < 2e-16 ***
## month.October     34.367      1.396  24.622  < 2e-16 ***
## month.November    20.042      1.396  14.359  < 2e-16 ***
## month.December     7.033      1.396   5.039 1.51e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.419 on 132 degrees of freedom
## Multiple R-squared:  0.9712, Adjusted R-squared:  0.9688 
## F-statistic: 405.1 on 11 and 132 DF,  p-value: < 2.2e-16

Обратим внимание, что январь пропущен. И февраль интерпертируется как разность между январем и февралем, март - разность между мартом и январем и т.д.

Косинусоидальный тренд

Пусть модель

\(\mu_t=\beta cos(2\pi ft+\Phi)\)

здесь \(\beta\) - амплитуда, а \(f\) - частота, естественно \(1/f\) период. Наблюдения повторяются через каждые \(1/f\) момента времени. В модель параметры \(\beta\) и \(\Phi\) входят нелинейно, но после преобразования

\(\beta cos(2\pi ft+\Phi)= \beta_1cos(2\pi ft)+\beta_2sin(2\pi ft)\)

где \(\beta=\sqrt{\beta_1^2+\beta_2^2}\) \(\Phi= atan(-\beta_2/\beta_1)\) Уже можно проводить линейную регрессию просто \(cos(2\pi ft),sin(2\pi ft)\) играют роль известных предикторов Подгоним модель тренда

\(\mu_t =\beta_0+\beta_1 cos (2\pi ft) +\beta_2 sin(2\pi ft)\)

к данным по температуре

har.<-harmonic(tempdub,1)
model4<-lm(tempdub~har.)
summary(model4)
## 
## Call:
## lm(formula = tempdub ~ har.)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1580  -2.2756  -0.1457   2.3754  11.2671 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      46.2660     0.3088 149.816  < 2e-16 ***
## har.cos(2*pi*t) -26.7079     0.4367 -61.154  < 2e-16 ***
## har.sin(2*pi*t)  -2.1697     0.4367  -4.968 1.93e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.706 on 141 degrees of freedom
## Multiple R-squared:  0.9639, Adjusted R-squared:  0.9634 
## F-statistic:  1882 on 2 and 141 DF,  p-value: < 2.2e-16

Здесь \(t\) измеряется в годах и частота наблюдений 1 год. На графике косинусоидальная модель и данные по температуре выглядят следующим образом

 plot(ts(fitted(model4),freq=12,start=c(1964,1)),
ylab="Temperature",type='l',col = "blue",lwd = 3,
ylim=range(c(fitted(model4),tempdub)))
points(tempdub,col="red",lwd=3)

Оценка дисперсии регрессионных оценок

Итак мы рассмотрели модель линейной регрессии вида

\(y_t=\mu_t+x_t\)

где стационарный процесс \(x_t\) имеет \(E[x_t]=0\),а \(\mu_t\) детерминированный тренд одного из рассмотренных выше типов: линейный, сезонный, косинусоидальный.

Начнем с простейшей сезонной модели. Оценки параметров \(\beta_1,\beta_2,...,\beta_{12}\) по сути есть средние за все январи, феврали,…, декабри. Если количество полных годов мы обозначим за \(N\) то среднее за каждый \(j\)-ый месяц года \(j = 1,..,12\) и есть оценка \(\widehat{\beta}_j\)

\(\widehat{\beta}_j=\frac{1}{N}\sum_{i=0}^{N-1}y_{j+12i}\)

Чтобы вычислить дисперсию оценки вспомним выражение (1). Заменяя в нем \(n\) на \(N\), а \(\rho_k\) на \(\rho_{12}\) получим

\(D[\widehat{\beta}_j]= \frac{c_0}{N}[1+2\sum_{k=1}^{N-1}(1-\frac{k}{N})\rho_{12}]\) для \(j=1,..,12\)

Если предположить, что \(x_t\)- белый шум, а следовательно все \(\rho_k=0\) то выражение выше сократится до

\(D[\widehat{\beta}_j]= \frac{c_0}{N}\) для \(j=1,..,12\)

Более того оно останется таким же, если несколько \(\rho_k\ne0\) , но \(\rho_{12}=0\)

Теперь посмотрим косинусоидальный тренд

Оценки параметров \(\beta_1\) и \(\beta_2\) амплитуд косинуса и синуса для всех частот вида \(m/n\) где \(m=1,2,...,n/2-1\)

\(\widehat{\beta}_1=\frac{2}{n}\sum_{t=1}^{n}\frac{cos(2\pi mt)}{n}y_t\) , \(\widehat{\beta}_2=\frac{2}{n}\sum_{t=1}^{n}\frac{sin(2\pi mt)}{n}y_t\)

Поскольку \(y_t\) в оценки входит линейно

\(D[\widehat{\beta}_1]=\frac{2c_0}{n}[1+\frac{4}{n}\sum_{s=2}^n\sum_{t=1}^{s-1}cos(\frac{2\pi mt}{n})cos(\frac{2\pi ms}{n})\rho_{s-t}]\)

где мы использовали факт,что \(\sum_{t=1}^n[cos(2\pi mt/n)]^2=n/2\)

Если \(x_t\)- белый шум, то все \(\rho_k=0\) и опять выражение для дисперсии сократится до \(D[\widehat{\beta}_j]= \frac{c_0}{N}\) Теперь рассмотрим модель с линейным трендом

\(\mu_t=\beta_0+\beta_1t: t=1,2,...,n\)

С оценкой параметра \(\beta_1\)

\(\widehat{\beta}_1=\frac{\sum_{t=1}^n(y_t-\overline y)(t-\overline t)}{\sum_{t=1}^n(t-\overline t)^2}\)

Здесь также оценка есть линейная комбинация \(y_t\), поэтому

\(D[\widehat{\beta}_1]=\frac{12c_0}{n(n^2-1)}[1+\frac{24}{n(n^2-1)}\sum_{s=2}^n\sum_{t=1}^{s-1}(t-t\overline)(s-t\overline)\rho_{s-t}]\)

Здесь использован тот факт, что \(\sum_{t=1}^n (t-t\overline )=\frac{n(n^2-1)}{12}\)

Интерперетация статистик вывода

Мы производили оценку регрессионных моделей при минимальных предположениях относительно шума \(x_t\). Однако вывод стандартных методов регресиионного анализа зависит от дополнительных предположений относительно распределения процесса \(x_t\). Стандартным предположением здесь является то, что \(x_t\) имеет нормальное распределение. Рассмотрим вывод результатов оценки параметров линейной регрессии для случайного блуждания

model1<- lm(x ~time(x))
summary(model1)
## 
## Call:
## lm(formula = x ~ time(x))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.041  -2.301   1.514   4.205   8.890 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.172910   0.789024  -6.556 4.65e-10 ***
## time(x)     -0.058242   0.006774  -8.598 2.40e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.572 on 199 degrees of freedom
## Multiple R-squared:  0.2709, Adjusted R-squared:  0.2672 
## F-statistic: 73.92 on 1 and 199 DF,  p-value: 2.396e-15

Процесс \(x_t\) имеет постоянную дисперсию. Оценивая дисперсию остатков, мы оцениваем дисперсию \(x_t\).

\(s^2 = \frac{1}{n-p}\sum_{t=1}^n(y_t-\widehat{\mu}_t)^2\)

здесь \(p\) количество оцениваемых параметров и \(n-p\) так называемое число степеней свободы для оценки стандартного отклонения \(s=\sqrt(s^2\)). Величина \(s\) постейшая мера качества оценки модели. Чем \(s\) меньше, тем лучше качество оценки. Другой мерой качества оценки есть величина \(R^2\) также называемая коэффициентом детерминации. Величина \(R^2\) интерпретируется как квадрат выборочной rкорреляции между наблюдаемыми данными \(y_t\) и оцененным трендом \(\widehat{\mu}_t\). Выраженный в процентах \(R^2 100\)% - сколько процентов изменчивости в наблюдениях объясняет подогнанный линейный тренд. Этот показатель очень удобен при сравнении моделей с разным числом параметров. Следует выбрать ту модель, где показатель \(R^2\) больше. Adjusted R-Squared интерпретируется также, но вычисляется по чуть измененной формуле. Стандартные отклонения оценок параметров помечены меткой Std.Error. С интерпретацией этих параметров следует быть аккуратным. Величины t-values это оценнный параметр модели, поделенный на его стандартное отклонение. В слечае нормальное распределенности \(x_t\) эти величины имеют\(t\)-распределение Стьюдента и используются для проверки нулевой гипотезы о том, что соответствующий параметр равен 0, а следовательно напрасно включен в модель. Подстановка значений t-values в фукцию t-распределения Стьюдента дает p-значения p-values. Что позволяет сделать выбор между отклонением или не отклонением нулевой гипотезы обоснованным для выбранного уровня значимости
F-statistics - статистика для проверки гипотезы об отсутствии связи между \(y_t\) и \(\mu_t+x_t\) при конкретном задании модели тренда \(\mu_t\). За ней следуют число степеней свободы F распределения вероятностей и Р-значение для проверки указанной выше гипотезы

Множественная линейная регрессия

Предположим, что _t представимо в виде

\(\mu_t=\beta_0+\beta_1 x_{t,1}+...\beta_{p}x_{t,p}\)

Где \(x_{t,k},k=1,...,p\) - известные постоянные. Обычно в качестве таких переменных берут величины , которые находятся под контролем экспериментатора и измеряются с пренебрежимо малой ошибкой, а \(\beta_0,\beta_1,...,\beta_{p}\) - неизвестные параметры, подлежащие оценке. Пусть значения \(x_{t,k}\) изменяются и исследователю доступны \(n\) их различных наборов и при этом наблюдается \(n\) значений \(Y_1₁,Y_2,...,Y_{n}\) переменной \(Y\), тогда

\(Y_{t}=\beta+\beta x_{t,1}+...\beta_{p}x_{t,p-1}+\epsilon_{t}, t=1,...,n\)

Или в матричной форме

\(Y=X\beta+\epsilon\)

Здесь мы предположили, что \(x_{t,0}=1, i=1,...n\)

Maтрица \(X\) размера \((n \times p+1)\) называется регрессионной матрицей. Переменные \(x_{t,k},k=1,...,p\) обычно называют независимыми переменными, регрессорами или предикторами, а переменнуя \(Y_t\) называют зависимой переменной или откликом.

Пример 1. Пусть \(x_{t,k}=t^{k}\). Получаем полиномиальную модель

\(Y_{t}=\beta_0+\beta_1 t+...\beta_{p}t^{p}+\epsilon_{t}, t=1,...n.\)

Пример 2. Пусть \(t\)- интерпретируется как время и принимает значания \(1, 2, ...n.\) \(x_{t,k}=x_{t-k}\) лаговая переменная, т.е. перееменная с временной задержкой в к единиц времени, тогда

\(Y_{t}=\beta_0+\beta_1x_{t-1}+...\beta_{p}x_{p}+\epsilon_t, t=1,.n\)

Последнюю модель называют моделью авторегрессии порядка \(p\). Общим в этих примерах является то, что они линейны по отношению к неизвестным параметрам \(\beta_0,\beta_1,...,\beta_{p}\) Поэтому модель называют \(линейной\) регрессионной моделью. Наиболее популярным методом получения оценок параметров \(\beta_0,\beta_1,...,\beta_{p}\). является так называемый метод наименьших квадратов. Метод заключается в минимизации суммы квадратов ошибок \(\epsilon′\epsilon=\sum_{i}\epsilon_{i}^2\) относительно вектора параметров \(\beta=[\beta_0,\beta_1₁,...,\beta_p]\)

Представим \(\epsilon′\epsilon\) в виде

\(\epsilon′\epsilon=(Y-X\beta)′(Y-X\beta)=Y′Y-2\beta′X′Y+\beta′X′X\beta\)

Продифференцируем \(\epsilon′\epsilon\) по \(\beta\). Пpиравняем полученную производную нулю получим

\(-2X′Y+2X′X\beta=0\)

или

\(X′X\beta=X′Y\)

Последние уравнения называют нормальными уравнениями. Решение этого уравнения вектор \(\tilde\beta\) называют оценкой наименьших квадратов. Простой проверкой нетрудно убедиться, что \(\tilde\beta\) доставляет минимум для \(\epsilon′\epsilon.\epsilon′\epsilon\) обозначают как RSS. Вектор \(\widetilde{Y}=X\tilde\beta=\tilde\beta_0+\tilde\beta_1x_{t,1}+...\tilde\beta_{p}x_{t,p}\) называют регрессией, а вектор \(e=Y-\widetilde Y\) называют остатками.

$e=Y-\widetilde{Y}=Y-X\tilde\beta$

Отметим, что при известных из линейной алгебры условиях на матрицу \(X′X\) оценки \(\tilde\beta\) и ошибки \(e\) будут единственны.

Свойства оценок наименьших квадратов

  1. Предположим, что вектор ошибок имеем нулевое математическое ожидание \(E[\epsilon]=0\). Тогда \(E[\tilde\beta]=(X′X)^{-1}X′E[Y]=(X′X)^{-1}X′X\beta=\beta\)

Таким образом \(\tilde\beta\) являются несмещеными оценками вектора \(\beta\).

  1. Предположим, что все \(\epsilon_{i} i=1,...,n\) некоррелированы и имеют одинаковую дисперсию, т.е. \(D[\epsilon]=\sigma^2I_{n}\) тогда

\(D[Y]=D[Y-X\tilde\beta]=D[\epsilon]=\sigma^2I_{n}\)

И

\(D[\tilde\beta]=D[(X′X)^{-1}X′Y]=(X′X)^{-1}X′D[Y]X(X′X)^{-1}=\sigma^2(X′X)^{-1}X′X(X′X)^{-1}=\sigma^2(X′X)^{-1}\)

Возникает вопрос почему именно \(\tilde\beta\) наиболее популярна при оценке параметров в линейной регрессионной модели. Ответ на этот вопрос дает следующая теорема

\(Теорема.\)

Пусть \(θ\) - оценка наименьших квадратов вектора \(θ=X \tilde\beta\). Тогда в классе всех других линейных несмещенных оценок линейной комбинации \(c′θ\), оценка \(c'θ\) является единственной оценкой с минимальной дисперсией.

До сих пор мы не делали никаких предположений о распределении. Если предоположить независимость и нормальность распределения ошибок \(\epsilon_t\sim N(0,\sigma^2)\) то Rao в 1974 году доказал, что \(c′\tilde\beta\) имеет минимальную дисперсию не только в классе всех линейных несмещенных оценок, но и в классе всех несмещенных оценок (эффективность). В этом случае также оценка максимального правдоподобия совпадает с оценкой наименьших квадратов. В случае не нормальности рапределения ошибок условия при которых оценка наименьших будет является ассимпототически эффективной найдены были Cox в 1968 году.

Дополнительно будем предполагать нормальность и независимость ошибок \(\epsilon_{i} \sim N(0,\sigma^2) или \epsilon∼N(0,I_{n}\sigma^2)\)

Теорема 2.Если \(Y\sim N(X\beta,I_{n}\sigma^2)\) и матрица \(X\) размера \((n⊗p)\) имеет ранг \(p\), тогда

  1. \(\tilde\beta \sim N(\beta,\sigma^2(X′X)^{-1})\)

  2. \((\tilde\beta-\beta)′X′X\tilde(\beta-\beta)/\sigma^2∼\chi^2\)

  3. \(\tilde\beta\) не зависит от \(s^2\) как случайные величины

  4. \(RSS/\sigma^2=(n-p)s^2/\sigma^2∼\chi_{n-p}²\)

Пример. Полиномиальная регрессия.

Готовим предиктор

q <- seq(from=0, to=20, by=0.1)

Вычисляем значения полинома.

y <- 500 + 0.4 * (q-10)^3

Генерируем нормальный шум, добавляем его к полиному и рисуем.

noise <- rnorm(length(q), mean=10, sd=80)
noisy.y <- y + noise
plot(q,noisy.y,col='deepskyblue4',xlab='q',main='Observed data')
lines(q,y,col='firebrick1',lwd=3)

Производим оценку модели и смотрим результаты

model <- lm(noisy.y ~ poly(q,3))
summary(model)
## 
## Call:
## lm(formula = noisy.y ~ poly(q, 3))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -253.277  -45.940    5.812   53.522  200.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   504.72       5.57  90.614   <2e-16 ***
## poly(q, 3)1  1819.86      78.97  23.045   <2e-16 ***
## poly(q, 3)2    53.87      78.97   0.682    0.496    
## poly(q, 3)3   904.11      78.97  11.449   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.97 on 197 degrees of freedom
## Multiple R-squared:  0.7708, Adjusted R-squared:  0.7673 
## F-statistic: 220.9 on 3 and 197 DF,  p-value: < 2.2e-16

Доверительные интервалы для оценок

confint(model, level=0.95)
##                 2.5 %    97.5 %
## (Intercept)  493.7358  515.7049
## poly(q, 3)1 1664.1238 1975.5894
## poly(q, 3)2 -101.8634  209.6022
## poly(q, 3)3  748.3733 1059.8389

Фрейм визуализации полиномиальной регрессии.

Пример. Множественная линейная регрессия

Читаем временные ряды.

Data<-read.csv("C:/tmp/Ruble.csv",header = TRUE)
head(Data)
##          X318x6 Rub.Usd brent  sp500 GAZPROM   Gold EUR.USD
## 1 1/11/16 11:00 75.9331 32.97 1915.6   13523 1105.9 1.09005
## 2 1/11/16 12:00 75.5411 33.29 1928.5   13547 1100.0 1.08833
## 3 1/11/16 13:00 75.3761 33.47 1921.6   13507 1099.2 1.08847
## 4 1/11/16 14:00 75.4561 33.20 1916.4   13451 1103.6 1.09006
## 5 1/11/16 15:00 75.3684 33.26 1921.0   13422 1102.6 1.08987
## 6 1/11/16 16:00 75.2620 33.25 1918.0   13436 1103.2 1.08983

Оцениваем модель c константой \[y_t=\beta_0+\beta_1 x_{t,1}+....+\beta_p x_{t,p}+\epsilon_t\]

regr <-  lm(Data$Rub.Usd ~ Data$brent + Data$sp500 + Data$GAZPROM+Data$Gold+Data$EUR.USD, data=Data)
summary(regr)
## 
## Call:
## lm(formula = Data$Rub.Usd ~ Data$brent + Data$sp500 + Data$GAZPROM + 
##     Data$Gold + Data$EUR.USD, data = Data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1042 -0.6984 -0.2149  0.5746  3.9030 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.833e+02  1.104e+01  16.611  < 2e-16 ***
## Data$brent   -5.191e-01  5.827e-02  -8.909  < 2e-16 ***
## Data$sp500   -1.151e-02  4.221e-03  -2.727  0.00676 ** 
## Data$GAZPROM -1.464e-03  4.223e-04  -3.467  0.00060 ***
## Data$Gold     5.427e-02  5.932e-03   9.149  < 2e-16 ***
## Data$EUR.USD -9.880e+01  1.062e+01  -9.300  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.106 on 312 degrees of freedom
##   (291 observations deleted due to missingness)
## Multiple R-squared:  0.6616, Adjusted R-squared:  0.6561 
## F-statistic:   122 on 5 and 312 DF,  p-value: < 2.2e-16
plot(Data$Rub.Usd,type = "b",col='deepskyblue4',xlab='Time',main='Rub/Usd')
lines(regr$fitted.values,col='firebrick1',lwd=3)

Остатки после удаления регрессии

plot(regr$residuals ,type = "b",col='blue4',xlab='Time',main='Rub/Usd.Residuals.')

Читаем временные ряды доходностей.

DataRates<-read.csv("C:/tmp/RubleRates.csv",header = TRUE)
head(DataRates)
##          X318x6      Rub.Usd        brent        sp500      GAZPROM
## 1 1/11/16 12:00 -0.005162439  0.009705793  0.006734183  0.001774754
## 2 1/11/16 13:00 -0.002184241  0.005407029 -0.003577910 -0.002952683
## 3 1/11/16 14:00  0.001061344 -0.008066926 -0.002706078 -0.004145998
## 4 1/11/16 15:00 -0.001162265  0.001807229  0.002400334 -0.002155974
## 5 1/11/16 16:00 -0.001411732 -0.000300661 -0.001561687  0.001043064
## 6 1/11/16 17:00  0.003057320  0.001203008  0.003076121  0.000893123
##           Gold       EUR.USD
## 1 -0.005335021 -0.0015779090
## 2 -0.000727273  0.0001286370
## 3  0.004002911  0.0014607660
## 4 -0.000906125 -0.0001743020
## 5  0.000544168 -0.0000367016
## 6 -0.002356780 -0.0015874040

Повторим регрессионный анализ для доходностей

regrRates <-  lm(DataRates$Rub.Usd ~ DataRates$brent + DataRates$sp500 + DataRates$GAZPROM+DataRates$Gold+DataRates$EUR.USD, data=DataRates)
summary(regrRates)
## 
## Call:
## lm(formula = DataRates$Rub.Usd ~ DataRates$brent + DataRates$sp500 + 
##     DataRates$GAZPROM + DataRates$Gold + DataRates$EUR.USD, data = DataRates)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.014016 -0.001919 -0.000009  0.001706  0.026740 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.0000792  0.0002117   0.374    0.709    
## DataRates$brent   -0.3006449  0.0247424 -12.151  < 2e-16 ***
## DataRates$sp500   -0.3121544  0.0729216  -4.281 2.48e-05 ***
## DataRates$GAZPROM  0.0509342  0.0532328   0.957    0.339    
## DataRates$Gold    -0.0660787  0.1003431  -0.659    0.511    
## DataRates$EUR.USD -0.0960630  0.1473420  -0.652    0.515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003747 on 311 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5358, Adjusted R-squared:  0.5283 
## F-statistic: 71.78 on 5 and 311 DF,  p-value: < 2.2e-16

Оцениваем модель без константы и без явно лишних предикторов \[y_t=\beta_1 x_{t,1}+....+\beta_q x_{t,q}+\epsilon_t\]

regrRates1 <-  lm(DataRates$Rub.Usd ~ DataRates$brent + DataRates$sp500  -1, data=DataRates)
 summary(regrRates1)
## 
## Call:
## lm(formula = DataRates$Rub.Usd ~ DataRates$brent + DataRates$sp500 - 
##     1, data = DataRates)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0139712 -0.0018486  0.0000472  0.0016979  0.0266948 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## DataRates$brent -0.29528    0.02312 -12.772   <2e-16 ***
## DataRates$sp500 -0.25138    0.05935  -4.235    3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003738 on 315 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5324, Adjusted R-squared:  0.5294 
## F-statistic: 179.3 on 2 and 315 DF,  p-value: < 2.2e-16
plot(DataRates$Rub.Usd,type = "b",col='deepskyblue4',xlab='Time',main='Rub/Usd.Rates.')
lines(regrRates1$fitted.values,col='firebrick1',lwd=3)

Остатки после удаления регрессии

plot(regrRates1$residuals ,type = "b",col='blue4',xlab='Time',main='Rub/Usd.Rates Residuals.')

Анализ остатков

После оценки регрессионной модели можно получить оценку процесса \(x_t\), а именно

\(\widehat{x}_t=y_t-\widehat{\mu}_t\)

которые будем назвать остатками. Если модель тренда была выбрана корректно, то поведение процесса остатков должно быть адекватно предположениям относительно процесса \(x_t\) перед оценкой модели. Например если предполагалось,что процесс \(x_t\) белый шум, то процесс остатков \(\widehat{x}_t\) тоже должен быть белым шумом. Если это не выполняется, это означает, что выбранная модель тренда не учитывает всей специфики динамики ряда. Начнем анализ остатков на примере данных по температуре, к которым оценивалась сезонная модель. Посмотрим на график нормированных остатков после удаления сезонной модели

data(tempdub)
plot(y=rstudent(model2),x=as.vector(time(tempdub)),xlab="Time",ylab="Standardized Residuals",type = "o", col = "blue",lwd=3)

model2
## 
## Call:
## lm(formula = tempdub ~ month. - 1)
## 
## Coefficients:
##   month.January   month.February      month.March      month.April  
##           16.61            20.65            32.48            46.52  
##       month.May       month.June       month.July     month.August  
##           58.09            67.50            71.72            69.33  
## month.September    month.October   month.November   month.December  
##           61.02            50.98            36.65            23.64

Нанесем символы месяцев на график остатков

plot(y=rstudent(model2),x=as.vector(time(tempdub)),xlab='Time',
ylab='Standardized Residuals',type='l',col = "blue",lwd = 3)
points(y=rstudent(model2),x=as.vector(time(tempdub)),
pch=as.vector(season(tempdub)))

Соберем остатки по месяцам

plot(y=rstudent(model2),x=as.vector(fitted(model2)),
xlab="Fitted Trend Values",
 ylab="Standardized Residuals",type = "n")
points(y=rstudent(model3),x=as.vector(fitted(model3)),
pch=as.vector(season(tempdub)))

О распределении остатков можно сделать предположения, если посмотреть на гистограмму остатков.

hist(rstudent(model3),xlab='Standardized Residuals', col = "blue")

Вид гистограммы наводит на мысль о нормальном распределении остатков. Нормальность остатков может быть проверена различными способами. Начнем с графического, так называемого нормального квантиль-квантиль (QQ) графика. Для нормально распределенных случайных величин QQ график должен выглядеть как прямая. Для остаков сезонной модели QQ график будет таким.

qqnorm(rstudent(model3),col="red")

Для примера смоделируем выборку с распределением Стьюдента, с похожей плотностью, и посмотрим на QQ график

О том, как выглядит плотность распределения Стьдента по сравнению с плтоностью нормального распределения напомнит следующий фрейм.

library("fGarch")
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:TSA':
## 
##     kurtosis, skewness
## Loading required package: timeSeries
## Loading required package: fBasics
## 
## Rmetrics Package fBasics
## Analysing Markets and calculating Basic Statistics
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.org
stud<- rstd(200,mean= 0,sd=1,nu=3)
hist(stud,xlab='t-Student simulates sample', col = "blue")

qqnorm(stud,col = "blue")

Прекрасный тест на проверку нормальности это тест Шапиро-Вилка (Shapiro-Wilk). О том, как устроен тест можно прочитать, например, здесь htps://en.wikipedia.org/wiki/Shapiro%E2%80%93Wilk_test

 shapiro.test((rstudent(model3)))
## 
##  Shapiro-Wilk normality test
## 
## data:  (rstudent(model3))
## W = 0.9929, p-value = 0.6954

Для примера проверка на нормальность выборки с распределением Стьюдента

shapiro.test(stud)
## 
##  Shapiro-Wilk normality test
## 
## data:  stud
## W = 0.9464, p-value = 8.624e-07

Не менее важно проверить не только нормальность (при исходном предположении модели), но и независимость (некооррелируемость) остатков Здесь нам поможет выборочная автокорреляционная функция, введенная в первой лекции

acf(rstudent(model3))

Хорошо видно, что все автокорреляции лежат внутри полосы \(\pm 2/\sqrt{n}\), что не отвергает гипотезу \(\rho_k=0\) для всех \(k=\pm1,\pm2,...\)

В качестве другого примера анализа остатков рассмотрим случайное блуждание и оценка линии регрессии.

model1<- lm(x ~time(x))
plot(x , col="blue" , type ="l", lwd=3)
abline(model1)

Здесь остаки после удаления тренда есть

plot(y=rstudent(model1),x=as.vector(time(x)),
ylab="Standardized Residuals",xlab="Time",type="o",col="blue",lwd=3 )

ACF остатков после удаления линейного тренда

acf(rstudent(model1))

ACF остатков заметно выходит за полосу \(\pm 2/\sqrt{n}\) таким образом не является белым шумом. Для процессов с такой ACF рассмотрим модели позднее в курсе. Напоследок проверим остатки на нормальность, построив QQ график остатков и выполним тест Шапиро-Вилка

qqnorm(rstudent(model1),col="red")

 shapiro.test((rstudent(model1)))
## 
##  Shapiro-Wilk normality test
## 
## data:  (rstudent(model1))
## W = 0.90268, p-value = 3.47e-10